//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007, 2008 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __octree_hpp__
#define __octree_hpp__

#include <string>
#include <map>

#include "math3d.hpp"
#include "utils.hpp"
#include "threading.hpp"
#include "dynamics.hpp"

class OctreeBuilder;

struct CubicScalar
{
  prec_t f;
  prec_t dfdx, dfdy, dfdz;
  prec_t d2fdxdy, d2fdxdz, d2fdydz;
  prec_t d3fdxdydz;

  void clear () { f = dfdx = dfdy = dfdz = d2fdxdy = d2fdxdz = d2fdydz = d3fdxdydz = 0; }
  void add (const CubicScalar& cv)
  {
    f += cv.f;
    dfdx += cv.dfdx;
    dfdy += cv.dfdy;
    dfdz += cv.dfdz;
    d2fdxdy += cv.d2fdxdy;
    d2fdxdz += cv.d2fdxdz;
    d2fdydz += cv.d2fdydz;
    d3fdxdydz += cv.d3fdxdydz;
  }

  void mul (prec_t v)
  {
    f *= v;
    dfdx *= v;
    dfdy *= v;
    dfdz *= v;
    d2fdxdy *= v;
    d2fdxdz *= v;
    d2fdydz *= v;
    d3fdxdydz *= v;
  }

  void serialize (prec_t *f, prec_t *dfdx, prec_t *dfdy, prec_t *dfdz,
      prec_t *d2fdxdy, prec_t *d2fdxdz, prec_t *d2fdydz, prec_t *d3fdxdydz)
  {
    f[0] = this->f;
    dfdx[0] = this->dfdx;
    dfdy[0] = this->dfdy;
    dfdz[0] = this->dfdz;
    d2fdxdy[0] = this->d2fdxdy;
    d2fdxdz[0] = this->d2fdxdz;
    d2fdydz[0] = this->d2fdydz;
    d3fdxdydz[0] = this->d3fdxdydz;
  }

  void deserialize (prec_t *f, prec_t *dfdx, prec_t *dfdy, prec_t *dfdz,
      prec_t *d2fdxdy, prec_t *d2fdxdz, prec_t *d2fdydz, prec_t *d3fdxdydz)
  {
    this->f = f[0];
    this->dfdx = dfdx[0];
    this->dfdy = dfdy[0];
    this->dfdz = dfdz[0];
    this->d2fdxdy = d2fdxdy[0];
    this->d2fdxdz = d2fdxdz[0];
    this->d2fdydz = d2fdydz[0];
    this->d3fdxdydz = d3fdxdydz[0];
  }
};

struct CubicVector
{
  vect3d f;
  vect3d dfdx, dfdy, dfdz;
  vect3d d2fdxdy, d2fdxdz, d2fdydz;
  vect3d d3fdxdydz;

  void clear ()
  {
    f.clear();
    dfdx.clear();
    dfdy.clear();
    dfdz.clear();
    d2fdxdy.clear();
    d2fdxdz.clear();
    d2fdydz.clear();
    d3fdxdydz.clear();
  }

  void add (const CubicVector& cv)
  {
    f += cv.f;
    dfdx += cv.dfdx;
    dfdy += cv.dfdy;
    dfdz += cv.dfdz;
    d2fdxdy += cv.d2fdxdy;
    d2fdxdz += cv.d2fdxdz;
    d2fdydz += cv.d2fdydz;
    d3fdxdydz += cv.d3fdxdydz;
  }

  void mul (prec_t v)
  {
    f *= v;
    dfdx *= v;
    dfdy *= v;
    dfdz *= v;
    d2fdxdy *= v;
    d2fdxdz *= v;
    d2fdydz *= v;
    d3fdxdydz *= v;
  }

  void serialize (prec_t *f, prec_t *dfdx, prec_t *dfdy, prec_t *dfdz,
      prec_t *d2fdxdy, prec_t *d2fdxdz, prec_t *d2fdydz, prec_t *d3fdxdydz)
  {
    f[0] = this->f.x;
    f[1] = this->f.y;
    f[2] = this->f.z;
    dfdx[0] = this->dfdx.x;
    dfdx[1] = this->dfdx.y;
    dfdx[2] = this->dfdx.z;
    dfdy[0] = this->dfdy.x;
    dfdy[1] = this->dfdy.y;
    dfdy[2] = this->dfdy.z;
    dfdz[0] = this->dfdz.x;
    dfdz[1] = this->dfdz.y;
    dfdz[2] = this->dfdz.z;
    d2fdxdy[0] = this->d2fdxdy.x;
    d2fdxdy[1] = this->d2fdxdy.y;
    d2fdxdy[2] = this->d2fdxdy.z;
    d2fdxdz[0] = this->d2fdxdz.x;
    d2fdxdz[1] = this->d2fdxdz.y;
    d2fdxdz[2] = this->d2fdxdz.z;
    d2fdydz[0] = this->d2fdydz.x;
    d2fdydz[1] = this->d2fdydz.y;
    d2fdydz[2] = this->d2fdydz.z;
    d3fdxdydz[0] = this->d3fdxdydz.x;
    d3fdxdydz[1] = this->d3fdxdydz.y;
    d3fdxdydz[2] = this->d3fdxdydz.z;
  }

  void deserialize (prec_t *f, prec_t *dfdx, prec_t *dfdy, prec_t *dfdz,
      prec_t *d2fdxdy, prec_t *d2fdxdz, prec_t *d2fdydz, prec_t *d3fdxdydz)
  {
    this->f.x = f[0];
    this->f.y = f[1];
    this->f.z = f[2];
    this->dfdx.x = dfdx[0];
    this->dfdx.y = dfdx[1];
    this->dfdx.z = dfdx[2];
    this->dfdy.x = dfdy[0];
    this->dfdy.y = dfdy[1];
    this->dfdy.z = dfdy[2];
    this->dfdz.x = dfdz[0];
    this->dfdz.y = dfdz[1];
    this->dfdz.z = dfdz[2];
    this->d2fdxdy.x = d2fdxdy[0];
    this->d2fdxdy.y = d2fdxdy[1];
    this->d2fdxdy.z = d2fdxdy[2];
    this->d2fdxdz.x = d2fdxdz[0];
    this->d2fdxdz.y = d2fdxdz[1];
    this->d2fdxdz.z = d2fdxdz[2];
    this->d2fdydz.x = d2fdydz[0];
    this->d2fdydz.y = d2fdydz[1];
    this->d2fdydz.z = d2fdydz[2];
    this->d3fdxdydz.x = d3fdxdydz[0];
    this->d3fdxdydz.y = d3fdxdydz[1];
    this->d3fdxdydz.z = d3fdxdydz[2];
  }
};

/* Octree tricubic interpolator. */
class Octree
{
public:
  Octree ();
  ~Octree ();

  operator bool () { return false; }

  bool getValues (const vect3d& pos, prec_t *scalars, vect3d *vectors);

  struct cell_t;
  struct subcells_t;

  struct interp_t
  {
    vect3d corner;
    float *fact;
  };

  struct subcells_t
  {
    vect3d middle;
    cell_t *sc[8];
  };

  struct cell_t
  {
    subcells_t *subcells;
    interp_t *interp;
    struct cell_t *next;
    size_t index;
  };

  bool save (const char *path);
  bool load (const char *path);
  bool save (const std::string& path) { return save (path.c_str()); }
  bool load (const std::string& path) { return load (path.c_str()); }

  friend class OctreeBuilder;

  void clear ();

  std::string getPovrayScript();

  struct level_t
  {
    vect3d bs; // box size
    vect3d hbs; // half box size
    vect3d rbs; // reciprocal box size
  };

  size_t size () const;

private:
  vect3d corner, box, hbox, rbox;
  size_t xc, yc, zc;
  typedef std::vector<level_t> levels_t;
  levels_t levels;

  size_t nScalars;
  size_t nVectors;
  size_t nValues;

  typedef std::vector<cell_t*> hash_t;
  hash_t hash;

  allocator<interp_t, 1024> a_interp;
  allocator<cell_t, 1024> a_cell;
  allocator<subcells_t, 1024> a_subcell;
  allocator<float, 4096> a_fact;

  typedef std::map<void *, size_t> cmap_t;

  void addPovrayScript(std::string& str, cell_t *cell, size_t level);
  void save (FPWriter& w, cmap_t& cmap, cell_t *cell);
  void get_stats (cell_t *cell, size_t& nCells, size_t& nInterp, size_t& nSubcell,
      cmap_t& cmap);
};

/* Octree builder for tricubic interpolation. */
class OctreeBuilder
{
public:
  OctreeBuilder (const vect3d& c1 = vect3d(-1, -1, -1), const vect3d& c2 = vect3d(1, 1, 1),
      size_t nScalars = 1, size_t nVectors = 0, size_t limitMemoryMB = 100, prec_t smallestResolution = 0.0);
  virtual ~OctreeBuilder ();

  void addLimit (const vect3d& n);
  void setLimitRadius (prec_t limitRadius);

  void build (Octree& tree);
  void rebuild (Octree& tree);

  virtual void addValues (const vect3d& pos, CubicScalar *scalars, CubicVector *vectors) = 0;
  virtual bool inContact (const vect3d& pos) = 0;

private:
  vect3d c1, c2;
  size_t nScalars;
  size_t nVectors;
  size_t nValues;
  size_t limitMemoryMB;
  prec_t smallestResolution;
  prec_t limitRadius;

  struct values_str {
    values_str() : values(0) { }
    ~values_str () { delete [] values; }
    vect3d pos;
    prec_t *f, *dfdx, *dfdy, *dfdz, *d2fdxdy, *d2fdxdz, *d2fdydz, *d3fdxdydz;
    prec_t *values;
  };

  struct o8entry_t
  {
    Octree::cell_t *cell;
    vect3d c1, c2;
    values_str *vals[8];
  };

  typedef std::vector<vect3d> limits_t;
  limits_t limits;

  bool isInside(const vect3d& c1, const vect3d& c2);

  struct fieldtask : public Task
  {
    size_t part;
    OctreeBuilder *builder;
    std::vector<values_str*> *values_map_calc;
    void execute () { builder->run_calc_values (*values_map_calc, part); }
  };

  struct celltask : public Task
  {
    size_t part;
    OctreeBuilder *builder;
    std::vector<o8entry_t> *cells;
    vect3d box;
    void execute () { builder->run_calc_cells (*cells, part, box); }
  };

  void run_calc_values (std::vector<values_str*>& values_map_calc, size_t part);
  void run_calc_cells (std::vector<o8entry_t>& cells, size_t part, vect3d box);

  typedef std::map<vect3d, values_str> values_map_t;
  values_map_t values_map;
  std::vector<values_str*> values_map_calc;
  std::vector<o8entry_t> cur;
};

#endif // !__octree_hpp__

